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We construct models of rotating stars using the perturbative approach introduced by J. Hartle 
in 1967, and a set of equations of state proposed to model hadronic interactions in the inner core of 
neutron stars. We integrate the equations of stellar structure to third order in the angular velocity 
and show, comparing our results to those obtained with fully non linear codes, to what extent third 
order corrections are needed to accurately reproduce the moment of inertia of a star which rotates 
at rates comparable to that of the fastest isolated pulsars. 
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PACS numbers: PACS numbers: 04.40.Dg, 97.60.Jd, 26.60+c 



I. INTRODUCTION 

In this paper we show that the structure of a rapidly 
rotating neutron star can be described to a high level of 
accuracy by using the perturbative approach introduced 
by Hartle in 1967 ij, subsequently developed to third 
order in the angular velocity f2 0. 

There are several reasons why the perturbative ap- 
proach may be preferred to the direct numerical inte- 
gration of Einstein's equations for stellar structure. One 
is that the angular behaviour of the solution is given ex- 
plicitely in terms of Legendre polynomials of first kind, 
whereas the radial behaviour can be found by integrating 
linear differential equations in the radial variable; thus, 
instead of solving the nonlinear, coupled, partial differ- 
ential equations of the exact theory one simply integrates 
a set of ordinary, linear differential equations. Moreover, 
these equations need to be integrated only inside the star, 
because the exterior solution can be found explicitely in 
terms of known polynomials. A further element of inter- 
est is that the perturbative approach allows to evaluate 
to what extent a physical quantity changes at each order 
in f2, and to check whether the n'^-order contribution 
is needed to estimate a choosen quantity or whether the 
first order estimates, often used to interpret astronomical 
observations, are enough. For these reasons, the pertur- 
bative approach expanded to 2nd order in rotation rate 
has been used by several authors to study properties of 
rotating neutron stars 

II II IE II II B (see also the 
review [lfj and refereces therein). 

On the other hand the perturbative approach has a 
drawback, since it fails when the angular velocity ap- 
proaches the mass-shedding limit, i.e. when the star ro- 
tates so fastly that the gravitational attraction is not 
sufficient to keep matter bound to the surface; numeri- 
cal integration of the exact equations shows that in these 
conditions a cusp forms on the equatorial plane of the star 
[Til Il2l IT^j l , which cannot be reproduced by the first few 
orders of the perturbative expansions (we shall discuss 
this point quantitatively in section IlIII) . However, the 
rotation rate of all known pulsars is much lower than the 
mass-shedding limit for most equations of state (EOS) 



considered in the literature. Thus, unless we are specifi- 
cally interested in the mass-shedding problem, we can use 
the perturbative approach to study the structural prop- 
erties of the observed neutron stars. The problem is to 
establish how many terms in the expansion we need to 
consider to describe, say, the fastest isolated pulsar ob- 
served so far, PSR 1937+21, whose period is T = 1.56 
ms, and this is the question we plan to answer in this 
paper. The equations of structure for uniformly rotating 
stars have been developed to second order in the angu- 
lar velocity in .1, ; 3j and subsequently extended to third 
order in Q. With respect to the derivation of Q we in- 
troduce two novelties: 

- the third order corrections involve two functions, w\ 
and 1^3 . W\ gives a contribution to the moment of iner- 
tia, W3 affects the mass-shedding velocity. Both w\ and 
u>3, contribute to the dragging of inertial frames. In Q 
the solution for w\ was given explicitely outside the star 
in terms of Legendre polynomials of second kind; that 
expression contained some errors that we correct. In ad- 
dition, we complete the third order solution giving the 
explicit expression of W3 outside the star. 

- In a procedure was developed to construct families 
of constant baryonic mass solutions with varying angular 
velocity, which was applied to study the behaviour of the 
moment of inertia along such sequencies, for small val- 
ues of the angular velocities. This procedure cannot be 
applied to fastly rotating stars and furthermore, it was 
shown to be unstable when the stellar mass approaches 
the maximum mass for any assigned equation of state. In 
this paper we develop, and apply, a different algorithm 
that allows to describe a sequence of constant baryonic 
mass stars at any value of the angular velocity (smaller 
than the mass-shedding limit), and which is stable also 
in the maximum mass limit. 

In section [n] we shall briefly explain what are the cor- 
rections to the metric and to the thermodynamical func- 
tions that are needed to describe the structure of a rotat- 
ing star to third order in Q, and in the appendix we shall 
summarize the equations they satisfy inside and outside 
the star. 

In order to test the perturbative scheme, we shall com- 
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pare the results we find solving these equations for three 
selected EOS, with those of ref. ^| an d 01 where the 
exact Einstein equations have been integrated for the 
same EOS. To hereafter, the adjective "exact" referred to 
Einstein's equations will mean that the full set of equa- 
tions is solved numerically and non-perturbatively. This 
comparison will be discussed in section ITTT1 and will allow 
us to assess the validity of the third order approach. 

We shall then integrate the third order equations to 
model a neutron star composed of an outer crust, an 
inner crust and a core, each shell being described by 
a non viscous fluid which obeys equations of state ap- 
propriate to describe different density regions. There is 
a general consensus on the matter composition of the 
crust, for which we use two well established EOS: the 
Baym-Pethick-Sutherland (BPS) EOS [3 for the outer 
crust (10 7 <C p & 4 • 10 11 g/cm 3 ), and the Pethick- 
Ravenhall-Lorenz (PRL) EOS [16| for the inner crust 
(4- 10 11 &p-&2- 10 14 g/cm 3 ). 

For the inner core, p > 2 ■ 10 14 g/cm 3 , we use recent 
EOS which model hadronic interactions in different ways 
leading to different composition and dynamics. The main 
assumptions underlying these EOS will be outlined in 
section Hv] 



For a given EOS there are two kind of equilibrium con- 
figurations, ususally indicated as "normal" and "supra- 
massive" . They both refer to rotating stars with constant 
baryonic mass and varying angular velocity. Each normal 
sequence has, as a limiting configuration, a non rotating, 
spherical star. The supramassive sequences do not pos- 
sess this limit, and all members of a sequence have a 
baryonic mass larger than the maximum mass of a non 
rotating, spherical star. Our study is confined to stars 
belonging to the normal sequence. 

The results of our calculations will be discussed in sec- 
tion |V| 



II. THE METRIC AND THE STRUCTURE OF 
A ROTATING STAR 



In this section we introduce all the quantities that de- 
scribe the star to third order in the angular velocity f2. 
Following P, 0, we write the metric as 



ds 2 



^ [1 + 2h(r, 9)] dt 2 
2m{r,ff) 



(1) 



[r - 2M(r)] 
r 2 [1 + 2k(r, 9)} {d9 2 + 



dr 2 

2a [d(/)-w(r,9)dt} 2 } 



sm 



Here e^ r \ e A W and M(r) are functions of r and describe 
the non rotating star solution of the TOV equations (see 
Appendix A). The functions h(r,6), m(r,9), k(r,9) and 



w(r, 9) are the perturbative corrections 

h(r,9) = h (r) + h 2 (r)P 2 (9)+O{n 4 ) 



m(r, 9) 
k(r, 9) 

w{r, 9) 



m Q {r) + m 2 (r)P 2 {9) + 0(f! 4 ) 
k 2 (r)P 2 (9) 

Mr)-fea(r)] P 2 (0) + O(O 4 ) 
u>(r) + iwi (r) 

1 dP 3 (9) 



(2) 
(3) 
(4) 

(5) 



w 3 (r)- 



sin 9 d9 



P 2 (9) and P 3 (9) are the I = 2 and / = 3 Legendre polyno- 
mials. As described in 0, 0, Q , the function uj is of order 
f2, whereas (ho,h 2 ,mo,m 2 ,k 2 ,v 2 ) are of order f2 2 , and 
(wijivs) are of order fl 3 . In a similar way, the displace- 
ment that an element of fluid, located at a given (r, 9) 
in the non rotating star, experiences because of rotation 
can be expanded in even powers of f2 



Z = Zo(r)+Ur)P2(9)+0(n 4 ), 



(6) 



where both £,o( r ) and £, 2 (r) are of order fl 2 . We shall 
assume that matter in the star is described by a perfect 
fluid with energy momentum tensor 



TH» = ( e + P)u^ L u v + Pg>- 



(7) 



Since the fluid element is displaced, its energy density 
and pressure change; in a reference frame that is mo- 
mentarily moving with the fluid, the pressure variation 
is 

5P{r,9) = [e(r) + P(r)] [<fp<,(r) + Sp 2 (r)P 2 (9)} (8) 
and it can be expressed in terms of £ as follows 

1 dP] 



Sp (r) 
5p 2 (r) 



Mr) 
-6(r) 



e + P dr 
1 dP 
e + P dr 



(9) 
(10) 



where e(r) and P(r) are the energy density and the pres- 
sure computed for the non rotating configuration. The 
corresponding energy density variation of the fluid ele- 
ment is 

6e(r,9) = ^(e + P)[6p (r)+6p 2 (r)P 2 (9)} . (11) 

This last equation has been derived on the assumption 
that the matter satisfies a barotropic EOS. 

In appendix A we shall write the equations satisfied by 
the corrections to the metric functions and to the ther- 
modynamical variables to the various order in f2, with the 
appropriate boundary conditions. It should be stressed 
that these equations have to be numerically integrated 
only for r < R, where R is the radius of the non rotating 
star. For r > R the solution can be found analytically. 
Hereafter, r — R and M(R) will indicate, respectively, 
the radius and the mass of the non rotating star, found 
by solving the TOV equations for an assigned EOS. 
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Once the various functions have been computed as 
shown in Appendix A, we can evaluate how the star 
changes in shape due to rotation, the moment of inertia, 
the mass-shedding limit and all quantities one is inter- 
ested in. Here we shall briefly summarize what are the 
relevant physical quantities we shall consider. 

• Stellar deformation 

The effect of rotation described by the metric JU 
on the shape of the star can be divided in two con- 
tributions: 

A) A spherical expansion which changes the radius 
of the star, and is described by the functions ho 
and mo. 

B) A quadrupole deformation, described by the 
functions h 2 , v 2 and m 2 . 

As a consequence of these two contributions, that 
are both of order fl 2 and that will be indicated 
with the labels A and B, respectively, the radius 
of the star, the eccentricity and the gravitational 
mass change as follows. 

A fluid element on the surface of the spherical non 
rotating star, i.e. at (r = R,6), will be displaced 
by an amount SR given by 

5R = SR A + 6R B = £ (i?) + &(R)P 2 (6) (12) 

where £o(R) and £ 2 (R) can be computed by eqs. 
(|A16|I and i|A26p as shown in appendix A. 

Once we know SR we can compute, for instance the 
eccentricity of the star 



E c = 



(equatorial radius) 2 /(polar radius) 2 



1 



1/2 



(13) 



The change in the gravitational mass is given by 

R 3 



SM a 



m (R) 



(14) 



Notice that 5M grav does not depend on quadrupo- 
lar perturbations (fi2, v 2 , m 2 ), nor on the perturba- 
tions that are first (u>) and third order (wi,ui3) in 
Q. The reason is that such quantity is invariant 
under rotations of the system, and it is even for 
parity transformations (which change £1 — > — f2). 

• Dragging of inertial frames 

The function u>(r, 9) in the metric Q is responsi- 
ble for the dragging of inertial frames. It affects 
two important physical quantities, i.e. the mass- 
shedding limit (section III A() and the moment of 
inertia (section III C|) . 



A. The mass-shedding limit 

The mass-shedding velocity can be found using a pro- 
cedure developed by 0. Given a star with assigned 



baryonic mass which rotates at a given fi, let us consider 
an element of fluid which belongs to the star and is lo- 
cated on the surface at the equator; from the metric Q 
it is easy to see that it moves with velocity 



ybouna 

where 



(r,e)-H{r.D) 



9. 



Lu(r) + wi(r) - -u>3(r) 



e ^(r,0) = r 2[i + 2k{r,8)]sm 2 6. 



(15) 
(16) 



In general, the fluid element will not follow a geodesic of 
metric 

Using the geodesic equations we can also compute the 
velocity of a particle which moves on a circular orbit just 
outside the equator, in the co-rotating direction 



yfree _ 




(w' + w' l — §U> 3 ) 

(U)' + W[ — IW3) ±-£ 

^ e 2 



(17) 



2>| 1/2 



where a prime indicates differentiation with respect to 
r. In general, for an assigned Q smaller than the mass- 
shedding limit the two velocities V bound and V^ ree are 
different, and in particular 

ybound ^ yfree 

When fl increases, the two velocities converge to the same 
limit; when this limit is reached the fluid element on the 
surface will not be bound anymore and the star will start 
loosing matter from the equator. The value of fl for 
which v bound — V^ ree is the the mass-shedding limit 
and will be indicated as Sl ms . It should be stressed that 
ybound anc j yfree are com p U ted at the equatorial radius 
of the rotating configuration, i.e. 



R rot = R + Z (R)-h 2 (R) 



(18) 



In order to apply this procedure, we need to construct 
sequences of stellar models with constant baryonic mass 
and variable angular velocity. 



B. The algorithm to construct a sequence of 
constant baryonic mass solutions 

The baryonic mass of the star is defined as 



L 



M bar = I ^J~g~ «* e d 3 x, (19) 



where g is the determinant of the 4-dimensional met- 
ric, d 3 x is the volume element on a t = constant hyper- 
surface, and eo(^) = rnNnbarir) is the rest mass-energy 
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density, not to be confused with the energy density e(r). 
Since all known interactions conserve the baryon number, 
eo is a conserved quantity, i.e. 



(eoO ;M = 0. 
The expansion of M bar in powers of f2 is 



where 



M bar = + 5M bar + 0(n 4 ) 



-1/2 



(20) 



(21) 



<°>=W "drfl-^] r 2 e (r) (22) 



and 



SM b , 



R 

47T / r 2 dr ( 1 



'0 
m (r) 



2M(r) 



-1/2 



r-2M{r) 3 
+^(e + P)^ (r)} . 



+ \r 2 [n- oj{r)fe- u{r) 



(23) 
eo(r) 



In the following we shall briefly compare the procedure 
developed by Hartle in 0] to construct constant baryonic 
mass sequences and the procedure we use. 



Hartle 's 



1. Given a central energy density e[r = 0) and an as- 
signed EOS, the TOV equations are solved to find 
the non rotating configuration; the model is iden- 
tified by a baryonic mass M bar = M. 

2. For an assigned value of the angular velocity f2 the 
equations of stellar structure are solved to order 
O 2 , imposing that the correction to the pressure, 
fipo( r = 0) is different from zero. The value of 
Spo(r = 0) is then changed until the sought bary- 
onic mass is reached. 



2. Our procedure 

Our procedure to generate sequences with constant 
baryonic mass is the following. 

1. Same as point 1 of Hartle's procedure. 

2. For an assigned value of f2 we solve the equations 
of stellar structure up to order il 3 , for the different 
values of e(r = 0), by imposing that 5po(r = 0) = 0. 
Among these models we select the one with bary- 
onic mass M. 

If, in addition, we want to compute the mass-shedding 
velocity, having computed all metric functions for a given 
n we evaluate V bound and V free . If V bound = V free we 



have reached the mass-shedding limit and we stop the 
procedure. If V bound < V^ ree we choose a higher value 
of ft and go to point 2. This procedure can be applied to 
rapidly rotating stars and to models of stars close to the 
maximum mass. 

Hartle's procedure presents a problem when applied 
to rapidly rotating stars; indeed if the rotation rate is 
sufficiently high the change in the central energy density 
Se(r — 0) can be so large that treating this change as 
a perturbation is inappropriate and produces a large er- 
ror in all quantities. With our approach, instead, this 
change is treated as a background quantity (we vary 
e(r = 0)), and the results are much more accurate (as can 
be seen by the comparison with the "exact" integrations 
of 0]>0I discussed in section ITTTjl . A second problem 
arises when the mass of the star is close to the maxi- 
mum mass (where we mean the maximum mass along a 
constant Q sequence). Since M bar = M bar (Q, e c ), where 



e(r = 0) + Se(r = 0) 



it follows that 

on 



dAL 



bar 



M bar dfl 
Since near the maximum mass 
8M bar 



'dM, 



bar 



it follows that 



de 
de c 

on 



£2 



M ba 



0, 



oo. 



Consequently, since in Hartle's procedure e(r = 0) is kept 
constant, and Se(r — 0) is treated as a perturbation, the 
procedure fails. 



C. The moment of inertia of a rotating star 

The angular momentum of an axisymmetric matter 
distribution rotating about the symmetry axis can be ex- 
pressed in terms of a conserved four vector, the angular 
momentum density 

where £^ is the Killing vector associated to axial sym- 
metry. The conserved total angular momentum J tot can 
be found by integrating J M over any spacelike hypersur- 
face. Given the coordinate system we use, the natural 



choice is a t = const hypersurface, and since 



J' 



-g J° d 3 x 



-gT\ d 3 x 



(24) 

where, as in (|19|l . g is the determinant of the 4- 
dimensional metric, and d 3 x is the volume element on 
the chosen hypersurface. 
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Consequently, the moment of inertia of the matter dis- 
tribution is 



/ = £ = ! 
n 



^ -lt=const 



-gT\ d A x. 



(25) 



Using Einstein's equations, in (|24[1 and (|25|) we can re- 
place T 1 ^ with -^R l ^ and expand in powers SI™. By this 
procedure it is possible to show that only terms with n 
odd contribute to this expansion, and that the first terms 
are 



j tot = j + 8 j + o(n 5 ) 

where J is the first order contribution 
f 



J = 



G 



dm 
dr 



(26) 



(27) 



\r=R 



and SJ is the third order correction 



1 



a dm 
dr 



h 



to 



dw\ 
dr 

4(v 2 - h 2 ) 



2M 



-■— (h>2 H — — — 

5 \ r - 2M 



(28) 



In eqs. <|27H and l|28|l the functions za and j(r) are defined 
as follows 



v = S1/27T, and we choose three EOS labelled as follows: 
AU - which results from an approach based on nuclear 
many body theory |2ffl 

APR2 - (APRb in Tj]) which takes into account recent 
results on nuclear many body theory pH |22| 
L - obtained from relativistic mean field theory I23J . 
AU and L have been used both in ^| and in |14| : the 
results of the two papers agree to better than 1% for 
both EOS, except that for the estimate of the moment 
of inertia for the EOS L, for which the relative error is, 
at most, of the order of 1.3 %. APR2 has been used in 
[3 only. Since we interpolate the EOS tables with the 
same linear interpolation as in 

[3 (in [l| a 5th-order 
polynomial has been used), in the following tables our 
results for AU and L will be compared with those of • 
For each EOS we consider two stellar sequences with 
constant baryonic mass and varying angular velocity: 
Sequence A, such that the gravitational mass of the non 
rotating configuration is M = 1.4 M©. 
Sequence B, that corresponds to the maximum mass. 
The corresponding values of the baryonic mass are given 
m tabled 



TABLE I: The baryonic masses of the stellar models used to 
compare the results of the perturbative approach with those 
found by integrating the exact Einstein equations are given 
for the three considered EOS. The data in column 2 and 3 
correspond to a non rotating star with gravitational mass 



vj{r) = Q — w(r), 



(29) 



1 



2M(r) 



Thus, the corrections to the moment of inertia are 

I = I (a) +51 + 0(n 4 ), (30) 

(31) 



r(o) - L 

si ' 



61 =" 

n 



It should be noted that the angular momentum defined 
in H24f> coincides with that defined using the asymptotic 
behaviour of the metric as in 

[3 (Ch. 19). In appendix 
A we shall show how to find 5J from the asymptotic 
behaviour of the metric function w\. 



III. A COMPARISON WITH 
NON-PERTURBATIVE RESULTS 

In order to verify to what extent the perturbative 
scheme correctly describes the physical properties of a 
rapidly rotating neutron star, we compare the results 
of the integration of the equations of stellar structure 
expanded to order SI 3 with those of ref. and [3 - 
where the fully non linear Einstein equations have been 
integrated. 

To perform the comparison we study stellar sequences 
with constant baryonic mass and variable spin frequency 





M bar /M e 


AU 


1.58 


2.64 


APR2 


1.55 


2.69 


L 


1.52 


3.16 



The values of v are chosen as in 

[3 and [3- 

The results of the comparison are summarized in Ta- 
bles^] (sequence A) and lllll fseauence B), where we show 
the stellar parameters for different values of the spin fre- 
quency. The lines labelled CST and BS refer, respec- 
tively, to the data found by integratin g th e exact equa- 
tions of stellar structure in [13j and in [LJj; BFGM refer 
to our results found by integrating the equations of stel- 
lar structure perturbed to order O 3 for the same models 
and EOS. The spin frequency v is given in column 2, the 
central density e c in column 3, the moment of inertia / 
in column 4, the gravitational mass M grav in column 5, 
the equatorial radius R eq in column 6, the ratio between 
the spin frequency and the mass-shedding frequency (as 
computed in [13J and |14j ) in column 7. For each value 
of v and for each quantity, we give the relative error A. 

We stress that, as explained in Section II, the quanti- 
ties e c , M grav and R eq are expanded in powers SI™ with n 
even, therefore, as far as our calculations are concerned, 
they are computed up to order SI 2 . Conversely, the an- 
gular momentum J tot is expanded in odd powers of S7, 
and we compute it to order SI 3 . The moment of inertia 
/ (see eq. I25fl is the angular momentum divided by SI, 
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TABLE II: Stellar parameters for the sequence A (see text). 
Data are computed along a sequence of stellar models with 
constant baryonic mass corresponding to a non-rotating con- 
figuration of mass M = 1.4 Mq and varying spin frequency 
v — £l/2n (given in KHz in column 2). e c is the central den- 
sity in units of e* = 10 15 g/cm 3 ; the moment of inertia I is in 
units of 7» = W 45 gem 2 , the equatorial radius R eq is in km. 
For each value of v and for each quantity, we give the rela- 
tive error A. M grav is the gravitational mass in solar mass 
units. In the last column the ratio of the spin frequency and 
the mass-shedding frequency (as computed in jl3| and [l4|') 
is given. 



EOS AU (M bar = 1.58 M Q ) 





v (KHz) 




I/I* 


M grav /M Q 


Req (Km) 


vjv ma 


CST 





1.206 


1.160 


1.400 


10.40 





BFGM 





1.204 


1.154 


1.395 


10.39 




A 




0.2% 


0.5% 


0.4% 


0.1% 




CST 


0.476 


1.192 


1.191 


1.403 


10.59 


0.379 


BFGM 


0.476 


1.190 


1.182 


1.398 


10.56 




A 




0.2% 


0.8% 


0.4% 


0.3% 




CST 


0.896 


1.148 


1.292 


1.412 


11.22 


0.713 


BFGM 


0.896 


1.155 


1.254 


1.405 


11.02 




A 




0.6% 


2.9% 


0.5% 


1.8% 




CST 


1.082 


1.111 


1.386 


1.412 


11.87 


0.861 


BFGM 


1.082 


1.132 


1.298 


1.410 


11.34 




A 




1.9% 


6.3% 


0.1% 


4.5% 




CST 


1.257 


1.048 


1.566 


1.432 


14.44 


1.000 


BFGM 


1.257 


1.108 


1.344 


1.414 


11.73 




A 




5.7% 


14.2% 


1.3% 


18.8% 




EOS APR2 (M bar = 1.55 M Q ) 


BS 





0.995 




1.403 


11.55 





BFGM 





0.988 


1.310 


1.389 


11.58 




A 




0.7% 




1.0% 


0.3% 




BS 


0.554 


0.970 


1.396 


1.409 


11.99 


0.532 


BFGM 


0.554 


0.964 


1.373 


1.394 


11.97 




A 




0.6% 


1.6% 


1.0% 


0.2% 




BS 


0.644 


0.960 


1.425 


1.411 


12.18 


0.619 


BFGM 


0.644 


0.956 


1.395 


1.396 


12.12 




A 




0.4% 


2.1% 


1.0% 


0.5% 




BS 


0.879 


0.920 


1.551 


1.418 


13.05 


0.844 


BFGM 


0.879 


0.929 


1.466 


1.401 


12.63 




A 




1.0% 


5.5% 


1.2% 


3.2% 




BS 


1.041 


0.870 


1.737 


1.428 


15.04 


1.000 


BFGM 


1.041 


0.905 


1.526 


1.405 


13.13 




A 




4.0% 


12.1% 


1.6% 


13.0% 




EOS L (M bar = 1-52 M ) 


CST 





0.433 


2.127 


1.400 


14.99 





BFGM 





0.440 


2.120 


1.402 


14.99 




A 




1.6% 


0.3% 


0.1% 


0.0% 




CST 


0.283 


0.427 


2.194 


1.402 


15.30 


0.395 


BFGM 


0.283 


0.433 


2.182 


1.404 


15.27 




A 




1.4% 


0.5% 


0.1% 


0.2% 




CST 


0.505 


0.411 


2.378 


1.407 


16.22 


0.705 


BFGM 


0.505 


0.419 


2.307 


1.408 


15.95 




A 




1.9% 


3.0% 


0.1% 


1.7% 




CST 


0.609 


0.398 


2.548 


1.412 


17.17 


0.851 


BFGM 


0.609 


0.410 


2.385 


1.410 


16.43 




A 




3.0% 


6.4% 


0.1% 


4.3% 




CST 


0.716 


0.377 


2.881 


1.419 


21.25 


1.000 


BFGM 


0.716 


0.400 


2.473 


1.414 


17.06 




A 




6.1% 


14.2% 


0.4% 


19.7% 





TABLE HI: The stellar parameters for the sequence B, cor- 
responding to the maximum mass (see text) are given as in 
table El 



EOS AU (M bar = 2.64 M Q ) 





V I, i \ 1 1 Z J 




TIT 


A/f 1 A/f 
ivlgrav / Mq 


Iteg (IVni ) 




CST 


o 


3.020 


1.982 


2.133 


9.41 


o 


NOI 


u 


3 090 
o.uzu 


1 Q(S4 


2.126 






A 

1 




0% 
u.u /o 


Q% 
u.y /o 


u.o /o 


n 9% 

u.z /o 




CST 


0.602 


2.547 


2.077 


2.141 


9.74 


0.357 


NOT 


u.ouz 


Z.ODD 


9 n^R 

Z.UOo 


z. loo 


Q 71 




A 
1 




7% 


Q% 
u.y /o 


u.o /o 


U.O /0 




CST 


1.174 


2.101 


2.281 


2.170 


10.40 


0.697 


NOI 


1.174 


2.159 


2.218 


2.158 


10.19 




A 




2.8% 


2.8 % 


0.6 % 


2.0 % 




CST 


1.481 


1.813 


2.525 


2.201 


11.20 


0.879 


NOI 


1.481 


1.968 


2.344 


2.178 


10.58 




A 




8.5% 


7.2% 


1.0% 


6.0% 




CST 


1.625 


1.637 


2.766 


2.227 


12.10 


0.964 


NOI 


1.625 


1.881 


2.416 


2.187 


10.82 




A 




15.0% 


12.7% 


1.8% 


10.6% 




CST 


1.684 


1.532 


2.974 


2.246 


13.66 


1.000 


NOI 


1.684 


1.844 


2.448 


2.192 


10.93 




A 




20.0% 


17.7% 


2.4% 


20.0 % 






EOS APR2 (M 


bar = 2.69 Mq) 




BS 


n 
u 


Z.OUU 




2.205 


in io 
1U. Iz 


o 


DF \ji\L 


U 


O 7/1 Q 

z. / 4o 


2.219 


2.202 


lU.Uo 




A 
i 




^ 7% 




0.1% 


Q% 




BS 


n finct 
u.ouy 


9 QOO 


2.352 


2.217 


i n a k 

1U.4D 


0.408 


DF vjlvl 


u.ouy 


O OOO 
z.zyz 


2.347 


2.212 


1U.44 




A 
Z_i 




U.o/o 


0.2% 


0.2% 


U.l/o 




BS 


1.081 


1.900 


2.593 


2.243 


11.21 


0.723 


BFGM 


1.081 


1.958 


2.532 


2.234 


10.99 




A 




3.0% 


2.4% 


0.4% 


2.0% 




BS 


1.188 


1.800 


2.687 


2.253 


11.50 


0.795 


BFGM 


1.188 


1.884 


2.589 


2.241 


11.16 




A 




4.6% 


3.6% 


0.5% 


3.0% 




BS 


1.284 


1.700 


2.799 


2.263 


11.86 


0.859 


BFGM 


1.284 


1.819 


2.645 


2.248 


11.33 




A 




12.0% 


5.5% 


0.7% 


4.5% 




BS 


1.494 


1.400 


3.337 


2.307 


14.17 


1.000 


BFGM 


1.494 


1.679 


2.788 


2.264 


11.79 




A 




20.0% 


18.8% 


2.0% 


17.0% 








EOS L 


(M ter = 3.16Af Q 






CST 





1.470 


4.676 


2.700 


13.70 


o 


DF \jl\L 





1.486 


4.530 


2.662 


13.63 




A 
i 




1.6% 


3.0% 


1.4% 


0.5% 




CST 


0.352 


1.201 


4.959 


2.706 


14.20 


0.341 


DF \jl\L 


0.352 


1.220 


4.786 


2.668 


14.09 




A 
i 




1.6% 


3.5% 


1.4% 


0.8% 




L/b 1 


0.714 


0.955 


5.554 


2.733 


15.24 


U.oyz 


BFGM 


0.714 


0.990 


5.211 


2.689 


14.84 




A 




3.6% 


6.0 % 


1.6% 


2.6% 




CST 


0.909 


0.802 


6.292 


2.765 


16.57 


0.880 


BFGM 


0.909 


0.896 


5.545 


2.707 


15.46 




A 




8.4% 


11.8% 


2.1% 


6.7% 




CST 


0.997 


0.710 


7.024 


2.792 


18.05 


0.966 


BFGM 


0.997 


0.852 


5.724 


2.715 


15.82 




A 




20.0% 


18.5% 


2.8% 


12.4% 




CST 


1.032 


0.655 


7.659 


2.813 


20.66 


1.000 


BFGM 


1.032 


0.835 


5.800 


2.720 


15.98 




A 




27.0% 


24% 


3.3% 


23.0% 
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therefore in this paper it is the sum of two contributions: 
/(°), of order ft coming from J computed at order O, 
and SI, of order fl 2 coming from SJ computed at order 
Q ; in this sense we refer to J as to a quantity computed 
up to 3rd-order in Q. 

It should be noted that a difference of the order of 
a percent between quantities computed by integrating 
the exact or the perturbed equations can be attributed 
essentially to interpolation of the EOS that are given in 
a tabulated form. The discrepancies are indeed of that 
order when v = 0, with the only exception of e c of the 
maximum mass star for APR2. The larger discrepancy 
found in this case (A = 5.7%) can be explained noting 
that the APR2 model yields a M(e c ) curve featuring an 
extended plateau around M = M max . For example, a 
7% change of the central density, from 2.65 x 10 15 g/cm 3 
to 2.85 x 10 15 g/cm 3 , only leads to a variation of the fifth 
digit in the value of the mass. 

From Table [H] we see that for the stars of sequence 
A a critical quantity in the comparison bewteen the ex- 
act and perturbative approach is the equatorial radius of 
the star; indeed, while the relative error A = (R eX act — 
Rperturb^) I Rexact IS smaller than 4.5% for v 0.85 v ms , 
it becomes of the order of 19 — 20% at mass-shedding. 
The reason of this discrepancy is that when the spin fre- 
quency approaches the mass-shedding limit the numeri- 
cal integration of the exact equations shows that a cusp 
forms on the e quat orial plane, from which the star starts 
loosing matter . This cusp is hardly reproduced by a 
perturbative expansion, and in any event it cannot be re- 
produced by truncating the expansion to order SI 3 . The 
physical quantity which is more sensitive to the error 
on the radius at mass-shedding is the moment of iner- 
tia (~ 14%). From Table HTl it also emerges that the 
agreement between exact calculations and 3rd-order per- 
turbative approach is better for the EOS APR2, which 
is stiff er than AU and L. 

For the maximum mass sequences B, Table ITTT1 shows 
that the discrepancies are larger; the largest error is on 
the central density, and the reason is the following. In 
order to construct a constant baryonic mass sequence, 
for a fixed value of v the central density is changed un- 
til the required mass is reached; since approaching the 
maximum mass 



TABLE IV: The mass-shedding spin frequencies are given for 
sequences A and B (see text) and for the considered EOS. 



,. de c 
hm — 

M^M max dMb ar 



CO, 



(32) 



when M — M max the determination of e c becomes less 
accurate. 

In Table llVl we show the mass-shedding frequency com- 
puted by the perturbative and the exact approach for the 
sequences A and B: the values we find are systematically 
higher than those found with the exact codes. The dis- 
crepancy has to be attributed to the inaccuracy of the 
perturbative approach near mass-shedding as explained 
above. 

However, it should be stressed that the rotation rate 
of known pulsars is much lower than the mass-shedding 



f ms in KHz 


Sequence A 




EOS AU 


EOS L 


EOS APR2 


BFGM 


1.533 


0.876 


1.299 


CST 


1.257 


0.716 




BS 


1.237 


0.728 


1.041 


Acst 


22.0% 


22.3% 




Ass 


23.9% 


20.3% 


24.8% 


Sequence B 




EOS AU 


EOS L 


EOS APR2 


BFGM 


2.057 


1.277 


1.842 


CST 


1.685 


1.032 




BS 


1.682 


1.022 


1.494 


Acst 


22.1% 


23.7% 




A B s 


22.3% 


24.9% 


23.3% 



limit computed for most equation of state considered in 
the literature; for instance, the rotation rate of the fastest 
isolated pulsar observed so far is v ~ 641 Hz, and unless 
we are specifically interested in the mass-shedding prob- 
lem, the results given in Tables ITT1 and lllll show that the 
perturbative approach can be used to study the struc- 
tural properties of the observed rapidly rotating neutron 
stars. Indeed, for v ~ 641Hz we find that the results 
of the perturbative and exact calculations for the EOS 
AU and APR2 differ by at most ~ 2% for all computed 
quantities, even for the maximum mass. 

For the EOS L, differences are larger, mainly because 
v = 641 Hz is closer to the mass-shedding frequency than 
for the other EOS, and the perturbative approach be- 
comes less accurate for the reasons explained above. In 
addition, unlike the EOS AU and APR2, the EOS L is 
tabulated with very few points, expecially at high den- 
sities (only four points for e c > 10 15 g/cm 3 ), and this 
introduces further differences which depend on the inter- 
polation scheme. 

The results of Hartle's approach have been compared 
to the exact results by several authors , Q , HJ • In all 
these papers the perturbative approach is developed up 
to order J! 2 . 

|6J focuses on the calculation of the mass-shedding fre- 
quency for the maximum mass configuration Qu m , which 
is computed for a number of EOS to see whether the scal- 
ing law defined in 0] an d based on the results of the in- 
tegration of the exact equations, is satisfied. The authors 
found mass-shedding frequencies systematically larger by 
^10 — 15% than those obtained from the "exact" scaling 
law. 

In 0, Hartle's procedure was applied to construct 
models of stars rotating at mass-shedding for three values 
of the mass and for EOS different from those we consider. 
The results were compared with those of ^3] where the 
exact equations were integrated for the same EOS, find- 
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ing that Hartle's approach leads to results compatible 
to exact calculations down to rotational periods of the 
order of 0.5 ms. However, this result was subsequently 
disclaimed in , where it was shown that the equatorial 
radius computed by the exact codes of was sistem- 
atically underestimated, leading to larger values of Sl ms . 
For this reason the agreement between and [13 was 
good. The conclusion of ^3 is that Hartle's method to 
order 17 2 cannot be considered very accurate in configu- 
rations approaching the mass-shedding limit, as we also 
find. 

Finally, in |9( the second order Hartle's approach is ap- 
plied to construct sequencies of stellar models with con- 
stant Mbar and varying SI for the same five EOS used to 
integrate the exact codes of 0] • They used a procedure 
different from ours, matching the angular momentum and 
gravitational mass computed by the two approaches vary- 
ing v and e c . Therefore in their approach J (computed 
to order f2) and M grav are equal by construction. The 
comparison between the perturbative and the exact ap- 
proaches is focused on the quadrupole moment Q, which 
is the coefficient of the r~ 3 P2(cos9) term in Newtonian 
potential |3J. To order fi 2 , Q is 

where K is defined in cq. (|A22|) . They found a relative 
error which depends on the EOS and increases with the 
rotation rate. For rates of the order of the fastest pulsar 
they found errors of the order of 20%. Our calculations 
of Q for the EOS APR2, AU and L that are also used 
in [14J show a much better agreement with the exact 
results, with a relative error smaller than 2% for M grav = 
1.4 M(T>. 



IV. THE EOS 

As stated in Section 1, in this work we have employed 
different models of the EOS of neutron star matter in the 
region of supranuclear density, corresponding to the star 
core, while the description of the outer and inner crust 
provided by the EOS of Refs. 0] and 0], respectively, 
has been kept fixed. The four EOS considered in the core 
have been obtained within two approches, nuclear many 
body theory (NMBT) and relativistic mean field theory 
(RMFT), based on different dynamical assumptions and 
different formalisms. 

The NMBT approach rests on the premise that nuclear 
matter can be described as a collection of nonrelativis- 
tic nucleons, interacting through phenomenological two- 
and three-body forces. While suffering from the obvi- 
ous limitations inherent in its nonrelativistic nature, this 
approach is made very attractive by its feature of being 
strongly constrained by the data. The two-nucleon po- 
tential provides a quantitative account of the properties 
of the two-nucleon system, i.e. deuteron properties and 
~ 4000 accurately measured nucleon-nucleon scattering 



phase shifts [24| , while the inclusion of the three-nucleon 
potential allows one to reproduce the binding energies of 
the three-nucleon bound states [25|. as well the empir- 
ical equilibrium properties of symmetric nuclear matter 
|22| . Calculations of the energies of the ground and low- 
lying excited states of nuclei with mass number A < 10 
|26j |. whose results are in excellent agreement with the 
experimental values, have also shown that NMBT has a 
remarkable predictive power. 

Within RMFT, based on the formalism of relativistic 
quantum field theory, nucleons are described as Dirac 
fermions interacting through meson exchange. In its 
simplest implementation the dynamics is described by 
a scalar and a vector meson field |27| . This approach is 
very elegant, but leads to a set of equations of motion 
that turns out to be only tractable in the mean field ap- 
proximation, i.e. treating the meson fields as classical 
fields. The meson masses and coupling constant are then 
fixed fitting the equilibrium properties of nuclear matter. 

Both NMBT and RMFT can be generalized to take 
into account the possible appearance of strange baryons, 
produced in weak interaction processes that may become 
energetically favourable at high density. However, the 
inclusion of baryons other than protons and neutrons en- 
tails a large uncertainty, as little is known of their in- 
teractions and the available models are only loosely con- 
strained by few data. 

Of the four EOS we consider, those referred to as APR2 
(also used in Section II I II to compare our results with 
those of non-perturbative approaches), BSS1 and BSS2 
are based on NMBT, while the one labelled G240 has 
been obtained from RMFT. The APR2 and BSS1 EOS 
include nucleons only, and their differences are mainly 
ascribable to the approximate treatment of the three- 
nucleon interactions in the BSS1 model and to the pres- 
ence of relativistic boost corrections in the APR2 model. 
The BSS2 model is a generalization of BSS1 including 
the hyperons S~ and A . Finally, the G240 EOS includes 
the full baryon octet, consisting of protons, neutrons and 
S ^, A and E ± . 

The non rotating neutron star configurations corre- 
sponding to the APR2 and BSS1 EOS exhibit simi- 
lar mass-radius relations, with a maximum gravitational 
mass ^, 2 M Q . The appearance of strange baryons makes 
the EOS softer, thus leading to sizable changes in the 
mass-radius relation. This feature is clearly visible in 
the case of both the BSS2 and G240 EOS 28], and leads 
to maximum masses of ~ 1.22 M and ~ 1.55 M Q , re- 
spectively. 

Comparison with the experimental determinations of 
the masses of neutron stars located in binary radio pulsar 
systems, yielding a narrow distribution centered at 1.35 
±0.04 M Q j2^, suggests that the dynamics underlying 
the BSS2 model leads to too soft an EOS. On the other 
hand, all other EOS appear to be compatible with the 
mass of the most massive neutron star listed in Ref. |29l ] 
(~1.56M ). 
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V. RESULTS 

Using the above EOS to describe matter in the crust 
and in the stellar core, and using the procedure developed 
in section ITll the equations of stellar structure perturbed 
to order fl 3 have been integrated to construct sequencies 
of stellar configurations with constant baryonic mass and 
varying angular velocity. In our study we consider six 
values of Mb ar , plus the maximum mass model for each 
EOS. The gravitational mass is 



EOS APR2 



where 



M grav = M 



M = 4tt 



5M, 



grav 



0(fi 4 ), 



e(r) dr 



M/M e 



M bar /M Q 


APR2 


BBS1 


BBS2 


G240 


1.31 


1.19 


1.18 


1.19 


1.20 


1.50 


1.35 


1.34 




1.37 


1.56 


1.39 


1.38 




1.41 


1.68 


1.49 


1.47 




1.51 


2.10 


1.80 


1.79 






2.35 


2.00 


1.97 







TABLE VI: Maximum masses of the non rotating configura- 
tions. 



EOS 




M max /M Q 


APR2 


2.70 


2.20 


BBS1 


2.41 


2.01 


BBS2 


1.74 


1.22 


G240 


1.35 


1.55 



-P 1.6 



(34) 



(35) 



"M bar =2.35 



M bar =2.10 



M 



bar 



=1.68 



M har =1.56 
M ba b r a iL50 

M bar =1.31 



v=0 

v =0.57 (KHz) 
v =0.76 (KHz) 
v=1.19 (KHz) 



2.5e+15 



e c (g/cm- 



is the gravitational mass of the non rotating configuration 
and 5M grav is given by eq. (|14J) : M grav depends on the 
EOS and on the angular velocity. In Table |V] for each 
EOS we give the values of M for the non rotating model 
that corresponds to the selected M bar 



TABLE V: For each value of the baryonic mass in column 1, 
the corresponding values of the gravitational mass of the non 
rotating configuration are tabulated for the selected EOS. 



FIG. 1: (Color online) The gravitational mass for the EOS 
APR2 is plotted as a function of the central density for dif- 
ferent values of the spin frequency, evidentiating the curves 
with Mb ar = const. 



A. Mass-shedding limit 



In Fig. |2we plot the mass-shedding frequency, v ms , as 
a function of the gravitational mass for each EOS. The 
dashed horizontal line is the frequency of PSR 1937+21, 
the fastest isolated pulsar observed so far, the mass of 
which is presently unknown. 



N 









APR2 




BBS1 




BBS2 




G240 




PSR 1937+21 







The maximum masses (baryonic and gravitational) of 
the non rotating configurations are given in Table IVll for 
the selected EOS. 

An example of how the gravitational mass changes due 
to rotation is shown in Fig. ^ where we plot M grav (in 
solar mass units) as a function of the central density for 
different values of the spin frequency for the EOS APR2, 
evidentiating the sequencies with M bar = const. 



M grav /M o 

FIG. 2: (Color online) The mass-shedding frequency v ma is 
plotted for each EOS as a function of the gravitational mass. 
The horizontal line corresponds to the spin frequency of PSR 
1937+21, which is the fastest isolated pulsar observed so far 
and whose mass is presently unknown. 



We see that the values of v ms for G240 are sistem- 
atically lower than those of the other EOS, while those 
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for APR2 are sistematically higher. This behaviour can 
be understood from the data given in Table IViTl where, 
for each EOS, we tabulate the ratio of the gravitational 
mass of the star rotating at mass-shedding to the equato- 
rial radius, M™^ v /R™ s , for three values of the baryonic 
mass. For a given mass, the ordering of v ma in Fig. 
corresponds to the same ordering of M^ v /R^ s in Ta- 
ble lVIII showing that more compact stars admit a higher 
mass-shedding limit. 



TABLE VII: For each EOS and for three values of the bary- 
onic mass, we tabulate the ratio of the gravitational mass of 
the star rotating at mass-shedding to the equatorial radius, 
Mgravl ' RTq ■ In the last three columns, we give the ratio M/R 
for the non rotating star with the same baryonic mass. 



EOS BBS2 




non rotating 
mass shedding 





v = v ma 


v = 




grav I *^eq 


M/R 


M bar /M Q 


1.31 


1.56 


2.35 


1.31 


1.56 


2.35 


G240 


0.106 


0.127 




0.131 


0.163 




BBS1 


0.114 


0.136 


0.210 


0.143 


0.169 


0.267 


BBS2 


0.117 






0.154 






APR2 


0.124 


0.147 


0.223 


0.152 


0.179 


0.268 



EOS APR2 



jo 



In the last three column of Table fVlTl we also give the 
ratio M/R of the non rotating star with the same bary- 
onic mass, and it is interesting to see that rotation af- 
fects the compactness of a star in a way which depends 
on the EOS. For instance, for M^ar = 1.31M Q (for which 
all EOS admit a stable configuration), the non rotat- 
ing star with EOS BBS2 (M/R = 0.154) is more com- 
pact than that with EOS APR2 (M/R = 0.152), while 
near mass-shedding the compactness of the BBS2 star 
(M™ s av /R™ s = 0.117) is smaller than that of the APR2 
star (M™JR™ = 0-124). 

To further clarify this behaviour in Fig. [3] we plot the 
energy density inside the star as a function of the ra- 
dial distance for these two EOS in the non rotating case 
and at mass-shedding. It is useful to remind that the 
main difference between the EOS APR2 and BBS2 is that 
in the core of the BBS2 star there are hyperons. From 
the upper panel of Fig. [31 we see that when the BBS2 
star does not rotate hyperons are highly concentrated in 
the core, making the central density very large; however, 
since the EOS is soft, when the star rotates fastly matter 
is allowed to distribute throughout the star. Conversely, 
due to the stiffness of the APR2 EOS, in the correspond- 
ing star (Fig. flower panel) the matter distribution does 
not seem to be so affected by rotation. A similar analysis 
can be found in were different EOS were considered. 

By comparing our results to order f2 3 with those to 
order fl 2 we find that the third order contributions to 
the mass-shedding limit are actually negligible, smaller 
than 1%, independently of the EOS. 



non rotating 
mass shedding 



r (km) 

FIG. 3: (Color online) The energy density distribution inside 
the star is plotted as a function of the radial distance for the 
EOS APR2 and BBS2, and for M bar = 1.31. e is given in 
units of e, = 10 15 g/crr? . In both cases the continuous line 
refers to the non rotating configuration, the dotted line to the 
star rotating at mass-shedding. 



B. Estimates of the moment of inertia 

In Fig. 0] we plot the moment of inertia as a function 
of the spin frequency, for each EOS and for given values 
of the baryonic mass up to the maximum mass. Since we 
know, as discussed in section III, that the mass-shedding 
frequencies estimated by our perturbative approach are 
about 20% larger than those found by integrating the ex- 
act equations, we plot the data only for v < 80% v msi 
where v ms is the mass-shedding frequency we find. From 
Fig. 0] we see that the moment of inertia has a non- 
trivial behaviour as the baryonic mass increases. In- 
deed, along a 17 = constant sequence of stellar models, if 
Mi, ar <g; MJ™: X then / is an increasing function of Mi, ar , 



11 



EOS APR2 



(a) 



0.6 0.8 1 

v (KHz) 



PSR 1937+21 



M 



=1.31 M n 



bar - ' ' ""o ' 
M bar =1.50 M 
M bar =1.56 M . 
M bar =1.68 M 



M 



bar 



=2.10 M n 



M bar =2.35 M 
M bar =2.60 M 
M bar= M max 



EOS BBS1 



0.6 0.8 1 

v (KHz) 



(b) 



PSR 1937+21 
.31 M„ 



M bar =1. 



bar : 
bar : 
bar : 



.50 M 
=1.56 M 
=1.68 M n 



M bar =2.10M o 
=2.35 M n 



M 



bar 
M bar= M n 



EOS BBS2 



EOS G240 




0.2 
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FIG. 4: (Color online) The moment of inertia (in units of J» == 10 4j <? cm 2 ) is plotted as a function of the spin frequency, for 
the EOS APR2 (panel a), BBS1 (panel b), BBS2 (panel c) and G240 (panel d). For each EOS we consider a few values of the 
baryonic mass up to the maximum mass, and plot the data only for v < 80% v ma (see text). The vertical dashed line is the 
spin frequency of PSR 1937+21. 



while as Mb ar approaches the maximum mass / becomes 
a decreasing function of M\, ar . 

This behaviour can be understood as follows: / ~ 
MR 2 , and while in general keeping constant as the 
mass increases the radius decreases, approaching the 
maximum mass the mass remains nearly constant while 
the radius is continuosly decreasing. This is true for any 
EOS, but of course it is more evident for softer EOS like 
BBS2 and G240. 

In a recent paper |3(j it has been suggested that a 
measurement of the moment of inertia from pulsar tim- 
ing data will impose significant constraints on the nuclear 
EOS. The authors consider the newly discovered binary 
double pulsar PSR J037-3039, and construct models of 
the fastest pulsar (M = 1.337 M , v = 276.8 Hz) with 
different EOS. They compute the moment of inertia as- 



suming the star is rigidly rotating using the fully non lin- 
ear, numerical code developed in 0], and compare the 
results with those obtained using Hartle's perturbative 
approach developed at first order in the angular velocity, 
i.e. they compute the quantity given in eq. (|31l) which 
we call . They show that the agreement is very good, 
and this had to be expected since the spin frequency of 
the considered star is quite low. What we want to show 
now is that if we consider stars that rotate faster, as for 
example PSR 1937+21, the first order approach is inaccu- 
rate, and we need to consider the third order corrections. 
To this purpose, in Fig. [5] we plot for each EOS and for 
Mbar = 1-56 Mq the first order term (dashed line) 
and / = + SI (solid line) where 51 is the correction 
found by solving the perturbed equations to order f2 3 . 
For each EOS, the data are plotted for v < 80% v ms . 
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FIG. 5: (Color online) The moment of inertia, measured in 
units of h = 10 45 gem 2 , is plotted versus the spin frequency 
(in KHz) for stars with baryonic mass Mbar = 1.56Mq. The 
continuous lines refer to I computed by solving the stellar 
structure equations to third order in the angular velocity, 
whereas the horizontal dashed lines refer to the first order 
term The vertical dashed line is the spin frequency of 

PSR 1937+21. For each EOS, the data are plotted only for 
v < 80% Ums (see text). 



The vertical dashed line refers to the spin frequency of 
PSR 1937+21. It is clear from the figure that the third 
order correction becomes relevant at frequencies higher 
than ~ 0.5 KHz. For instance if v = 641 Hz the rela- 
tive error [51/(1^ + SI)] that one makes neglecting the 
third order corrections are 5% for APR2, 7% for BBS1 
and 13% for G240. 

It is interesting to compare 

/(o), I = l(°)+5I, where 51 
is computed using corrections of order O 3 , and the value 
of I one gets by integrating the exact equations of stellar 
structure. This is possible for the EOS APR2 and the 
results are shown in Table rVTTTl for M grav ~ 1.4 M© and 
v = 644 Hz, and for the maximum mass and v = 609 Hz. 

The table shows that including third order is impor- 
tant to adequately estimate the moment of inertia of the 
fastest pulsars. 



VI. CONCLUDING REMARKS 



TABLE VIII: We compare the moment of inertia computed 
for the EOS APR2 by integrating the exact equations of stel- 
lar structure (column 2) |14| . the perturbed equations to order 
Q (column 3), and to order Q, 3 (column 5). The relative error 
bewteen the exact and the approximated results are given in 
columns 4 and 6. I is in units of /» = 10 45 gem 2 . 





I exact 1 1* 


I^/h 


A 


(i^+Sl) /I, 


A 


M grav ~ 1.4 M Q 
v = 644 Hz 


1.425 


1.238 


13% 


1.395 


2.1% 




I exact 1 1* 


I (0) /h 


A 


(i^+Sl) /h 


A 


-Mgrav — -M m ax 

v = 609 Hz 


2.352 


2.281 


3% 


2.347 


0.2% 



approach fails to correctly reproduce the stellar proper- 
ties and the reasons are well understood. We confirm this 
behaviour and give quantitative results in this limit. 

However, for lower values of the angular velocity the 
situation is different. Taking as a reference the rotation 
rate of the fastest isolated pulsar observed so far, PSR 
1937+21, for which v ~ 641 Hz, we see that at these rates 
the perturbative approach allows to describe all stellar 
properties to an accuracy better than ~ 2% even for the 
maximum mass models, unless the EOS in the core is 
very soft, as in the case of the EOS L; indeed, for this 
EOS the mass-shedding velocity is low and close to the 
chosen rotation rate, so that the perturbative approach 
is inaccurate. 

The two quantities that are affected by third order cor- 
rections are the mass-shedding velocity and the moment 
of inertia. For the first, we find that the third order 
corrections are actually negligible, smaller than 1%, in- 
dependently of the EOS. 

Conversely, the moment of inertia is affected by terms 
of order O 3 in a significant way; for instance, for a star 
with M = 1.4 M© rotating at v ~ 641 Hz the third order 
correction is of the order of 51/(1^°' + SI) ~ 5% for the 
stiffer EOS we use (APR2), and as high as ~ 13% for the 
softest (G240). 

Thus, third order corrections have to be included to 
study the moment of inertia of rapidly rotating neutron 
stars, whereas they are irrelevant to estimate the mass- 
shedding limit. 



In this paper we have solved the equations of stellar 
structure developed to third order in the angular veloc- 
ity, using recent EOS which model hadronic interactions 
in different ways to describe the matter in the inner 
core. The stellar parameters we find have been com- 
pared, when available in the literature, with those found 
by solving the exact equations of stellar structure. The 
main conclusions we can draw from our study are the 
following. 

It is known that near mass-shedding the perturbative 



APPENDIX A 

In this appendix we write the equations that govern the 
radial part of the metric and thermodynamical functions 
that describe the structure of a rotating star to third 
order in the angular velocity O. For each function we 
specify the appropriate boundary conditions and we show 
how to construct the solution order by order. We give 
explicitely the analytic solution when available, following 
and completing the work of Hartle and collaborators [J, 
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0, 0j- In what follows we shall refer to the functions 
appearing in section ITT1 

We recall that the structure of a non rotating, spherical 
star is described by the TOV equations 



dM{r) 

dr 
dv{r) 

dr 
dP{r) 

dr 



4irr 2 e(r), 

} M(r) + 4nr 3 P(r) 
J r(r - 2M(r)) ' 
(P(r)+e(r)) 



(Al) 



e(r) and P(r) are the energy density and pressure distri- 
bution in the non rotating star. Since all functions we 
shall consider depend on r, we shall omit this dependency 
throughout. 

Hereafter, R and M{R) will indicate, respectively, the 
radius and the mass of the non rotating star, found by 
solving eqs. IjAlfl for an assigned EOS. 



1. The first order equations 

To order Q, the shape of the star remains spherical 
and there is only one function, lo, to determine, which is 
responsible for the dragging of inertial frames (cfr. the 
expansion in eq. EJ. By introducing the function vj = 
fl — lj, it is easy to show that it satisfies the following 
equation 



Id, dzu . 4 dj 
-{r J-j-) H 



dr 



dr 



• dr 



0, 



(A2) 



where j(r) = e u ' 2 yl — 2 y-. Equation \K2\ can be re- 
duced to two first order equations by defining two auxil- 
iary functions of r 



X = J w j 



r J 



.dm 

dr 



(A3) 



that satisfy 

r<R 



dr 
du 
dr 



47r,r 2 (e + P)x 



r 4 r-2M 
167rr 5 (e + P)x 



2M 



(A4) 



They have to be integrated from r 
following boundary conditions 



to R with the 



x(o)=j(oK 



u(0) = 0. 



(A5) 



These conditions follow from the behaviour of w near the 
origin 



w ~ vj c (\ + nj 2 r 2 + ....), 

where w c is a constant. Note that from eqs. 
6j) it fopllows that u goes to zero faster than r 4 . 



(A6) 



31 and 



For r > R, P and e vanish, 
solution of eqs. (|A4() is 



r > R 



X (r) - n 



2J 



M = M(R), and the 
u{r) = 6J, (A7) 



where J is a constant which represents the angular mo- 
mentum of the star to first order in Vt. The constants w c 
and J can be found by imposing that the interior and the 
exterior solutions match continuously at R, i.e. 



2J 

w c x(-R) = n - — , 



u(R) = 6 J. 



2. Second order equations 

The second order terms affect the structure of the star 
in the following way: ho and mo produce a spherical 
expansion, /12 and V2 a quadrupolc deformation. 

On the assumption that the EOS is a one parameter 
equation of state, Einstein's equations admit a first in- 
tegral (cfr. |l|) that allows to find two equations, one 
relating Hq and Sp , the second hi and <5p2- These two 
equations are 



|U = <5p + h - -e "r 2 w 2 
Q = 5p 2 + h 2 + \e- v r 2 m 2 , 



(A8) 
(A9) 



where /i is the correction of order fi 2 to the chemical 
potential /j, c 



tx c = M [i+M+o(r! 4 )] 



£ + V 



■ exp 





r d£ 


J 


' £ + V 



(A10) 



(where £ and V are the energy density and pressure in 
the rotating configuration) which can be shown to be 
constant with respect to r and 9 throughout the star. 



a. Spherical expansion 

Using eq. i|A8(l , the relevant equations for the spherical 
deformation can be shown to be 



4irr~ 



de_ 
dP 



M(e + P)} 



12r 4 



drrio 
dr 

dSpo u 2 
~dV = 12r 4 (r - 2M) 
mo(l + 8irr 2 P) 4n(e + P)r 2 Sp 



8irr 5 (e + P) X 2 
3(r - 2M) 



(r - 2M) 2 

2r 2 X 
"3(r- 2M) 



u 



(r 



r - 2M 

- 3M - 47rr 3 P)x 



r - 2M 



(All) 



where u and \ are defined in l|A3(l . These equations must 
be integrated inside the star from r = to R with the 
condition that both too and Spo vanish in r — 0. To start 
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the numerical integration it is useful to make use of the 
following expansion 



where the constants A and B are related by the following 
expression 



m (r) 



-in 
15 



e(0) + P(0)] 



Spo(r) ~ - x(0) r 2 , 



de_ 
dP 

0. 



X(0) 2 r 5 



(A12) 



For r > R u and \ reduce to (|A7|) : by integrating the 
first eq. (|A11|) with P and e vanishing, one finds 



m (r) 



J 



8M 



(A13) 



where (5M is a constant which represents the correction 
to the gravitational mass due to the spherical expansion 
of the star. This quantity can be determined by imposing 
the continuity of the solution at r = P, i.e. 



SM = m (R) 



p 3 ' 



(A14) 



and mo(-R) is found by integrating eqs. (|A11|) for r < R. 

Due to the sperical expansion, the pressure of each 
element of fluid changes by an amount Spo(r) found by 
integrating eqs. lAllf) : at the same time the element is 
radially displaced by an amount 



r + £o(r), 



(A15) 



and £,o(r) can be found in terms of Spo(r) using eq. J^J; 
known £o ( r ) we can compute the corresponding variation 
of the stellar radius (cfr. eq. IT^I 



6R A =£o(R) 



Spo(R)R(R-2M(R)) 
M(R) ' 



(A16) 



b. Quadrupole deformation 

We shall now solve the equations for h 2 (r) e v 2 (r), 
that are responsible for the quadrupole deformation of 
the star: 



dv 2 
dr 

dh 2 
dr 



dv , ,1 
■— h 2 + - 
dr r 



dv 
dr 
4v 2 



1 dv, 
2dr J 



87rr 5 (e + P)x 2 



3(r-2M) 



u 

6^ 



B + 2tt 



e(0)+P(0) 



A=^r[e(0) + P(0)] (j(0)n7 c ) 2 . 



(A19) 

Following [lj , we write the general solution of eqs. (|A17|) 
for r < R as 



h 2 (r) = h% + Chf 
v 2 (r) = uf + C vg. 



(A20) 



C is a constant to be determined (see below), (h 2 ,v 2 ) 
are a particular solution found by integrating eqs. (|A17|) 
with the initial conditions, say, A = 1 and B given by 
eq. (|A19|I : (h 2 ,v 2 ) are the solution of the homogeneous 
system which is found by putting \ = and u = in 
eqs. <|A17ll . numerically integrated with the boundary 
condition 







h¥(r) 



(A21) 



-2vr 



-e(0) + P(0) 



For r > R eqs. (|A17|l can be solved analytically in 
terms of the associated Legendre functions of the second 
kind, and the solution is 0,0 



r>R h 2 (r) = J 2 ( 



1 



J2 

v 2 (r) = j 



M(R)r 3 
K 



l-) + KQl(0 (A22) 



2M(R) 



where if is a constant and 
r 3 



QUO 
QUO 



L 2 



r[r-2M(R))]V 



£+F 3^-5^ 



;QU0 



a/2 



£-1' 

3^_2 

e 2 -i 



(A23) 



1. 



with 



M(P) 



(A24) 



2M v dr 
,dv 



M^ + p)-^) 



r(r - 2M) dr 

87rr 5 (e + P)x 2 
3(r - 2M) 



1 di/ 

2 dr r 



r" 
1 

~2M v dr 



-(-) 



(Al7)We can now determine the constants C and if by im- 
posing that the solution l|A20(l and the solution (|A22|I 
are equal at r = P; thus, the functions h 2 and v 2 are 
completely determined. 

We can now compute the quadrupole deformation of 
the star. Indeed, known h 2 from eq. l|A9fl we find 5p 2 : 



1 dv 

r H 

2 dr r 



1 



2M v dr 



dv 



^2^) = -/l2W 



1 



r w{rY 



(A25) 



Let us consider the solution for r < P first. 
It is easy to check that a regular solution of eqs. (|A17|I 
near the origin must behave as 



and from eq. (|1(J|) we find the second correction to the 
displacement 



r -> h 2 (r) ~ Ar 2 



U 2 ~ Pr 4 



(A18) 



6W = -feW/ 



1 dp 



P dr 



(A26) 
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Thus, an element of fluid located at a given r and a given 
9 in the non rotating star, when the star rotates moves 
to a new position identified by the same value of 6, and 
by a radial coordinate 



f = r + £o(r) + 6(r)P 2 (0) + 0(fi 4 ), 



(A27) 



where £o( r ) an d £,2( r ) are found as explained above. 

The remaining metric function to,2 can be determined 
by a suitable combination of Einstein's equations, which 
gives 



m 2 = (r - 2M) 



-hi 



8nr 5 {e + P)x 2 u 2 
3(r - 2M) + 6^ 



(A28) 



3. Third order equations 

The functions that describe the third order corrections 
are w\ and W3, and satisfy two second order differential 
equations with a similar structure. Since they do not 
couple, we shall solve them independently. 



a. The equations for wi 



C is a constant, (xf , wf ), are a particular solution found 
by integrating eqs. 1A32I) with the initial conditions 

X f (0) = 0, uf (0) = 0; 

(xf , u i) is the solution of eqs. (|A32|) in which D and D 2 
are set to zero, numerically integrated with the asymp- 
totic conditions 

r -> xf (r) ~ 1 + y [e(0) + P(0)] r 2 , (A34) 



«f(r) 



16vr 



s(0) + P(0)] r 5 



For r > P eqs. (IA32ft reduce to 



dxi _ «1 

dr r 4 ' 

dui u d 6 J 2 



(A35) 



Note that, since j(r > R) = 1, for r > R xi and wi 
coincide (cfr. eq. IA31H . 

At large distance from the source w\ must behave as 



/ v 26J n .l. 

Mr) = -r + o(-j) 



(A36) 



r 4 dr 
where 



1 d 4 .dwi(r) 4 dj 1 
^31 ( r J — ; ) H r"! = -Do - --D 2 , 



dr 



■ dr 



r 4 D = -u— 
dr 



to 



2toq 



■ + (" 



r - 2M 
1 



r - 2M v (dP/de) 
4fc 2 - /i 2 



r 4 £>2 u d 



5 5 dr 
2to 2 



- h 



l)S Po 



TO 2 



-4r 3 



'X 



1$ 
j dr 
2r 3 x 2 



1 



_r-2AT v (dP/de) 
By introducing the variables 



r-2M 



1)6P2 



3(r - 2M) 



(^) 



2r 3 x 2 



3(r - 2M) 



Xi = jwi 
eq. I|A29|) becomes 

dxi _ ui 47rr 2 (e + P)xi 



4 .du>i 

"1 = r J—r 
dr 



dr r 4 r - 2Af 

dtti _ 167rr 5 (e + P)xi 



(A29) 



(A30) 



I* 

j dr 



-i (A31) 



r 4 P - —D 2 . 
5 



dr r - 2M 

We write the general solution of this system for r 

Xi = xf + c xf 



(A32) 

< R as 
(A33) 



Ml = tif + C Zif . 



where J J is the correction to the angular momentum of 
the star; this suggests to write w±, and therefore xi> as 

2(5 J 

w 1 (r)=xi(r) = — +F(r), (A37) 

where F(r) must go to zero at infinity faster than r 3 . 
Consequently, when r — > 00 Hi — > — 6<5J. Thus, by 
direct integration of the second eq. (| A35|) . and using eqs. 
(|A7|) . © and ljA22|l . we find 



: [4fc 2 



6J^ 

r.4 



6(5 J 



84J 3 



24J 3 



24JiTQ 2 . 



r* " 5r 4 5M(P)r 3 

48 A/ (P) JPTQ^ 



5[r(r - 2M(P))] 1 /2 



- 6(5 J. 



(A38) 



Knowing u1 xt , we can now integrate the first of eqs. 
(|A35|I which, using eq. ljA37|) . becomes an equation for 
F(r), the solution of which is 



12J 3 



108r 4 ln(- 



4J 3 
5Afr 6 
r 



F(r) = - 
JK 

40M 3 r 4 |_~~"" ""V - 2M 
+33r 4 - 240r 3 Af + 336r 2 M 2 



) - 288r 3 Afln(- 



Y-2M 
256APY - 96M 4 



(A39) 
) 



-192rAf 3 ln( 



' ) + 12r 4 ln( ' 



r - 2M 



2 AT 



33JP 
40 A/ 3 ' 



where M = Af (P). This expression satisfies eq. (|A35|) . as 
it can be checked by direct substitution; it may be noted 
that it differs from the expression of F(r) given in ref. 
U which, conversely, does not satisfy ea. l|A35(l . 
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In conclusion 



xT 



2SJ 



F(r). 



(A40) 



At this point the situation is the following. We have in- 
tegrated eqs. (IA32|) for r < R, and we know the solution 
for i*i and \i up to the constant C (eq. IA33|I : we have 
found the solution for r > R for u\ fea lA38|l and for Xi 
(eq. IA40|) up to the constant 5 J. The two constants are 
found by imposing that the solutions for u\ and xi given 
by eqs. (|A33|I and by eqs. (|A38IA40|) match continuously 
at r = R. 

It should be noted that computing 5 J by this proce- 
dure or by eq. (|2*5|) is equivalent; indeed we find a relative 
difference smaller than 10 -10 . This is a further check of 
the correctness of the analytic expression (|A39|) for F(r) 
we derived. 



b. The equations for W3 

The solution for the remaining function W3 can be 
found along the same path. We start with the equation 



1 d ^ ri . dw z {r) | 4 dj ^ 
r 4 dr dr r dr r(r — 2M) 5 



10jW3 =iA,(A41) 



and assume r < R; we introduce the functions 



X3 = jw 3 

and find 

d-X3 _ u 3 47rr 2 (e + P)x 3 



u 3 = r^ (A42) 
dr 



dr 



r-2M 



(A43) 



dus = 167rr 5 ( £ + F)x3 Wxsr 3 r 4 ^ 
dr r-2M (r - 2M) 5 2 ' 

The general solution of eqs. I|A43(1 for r < R is 

X3 = Xa + C X3 (A44) 
U3 = U3 + C u 3 . 

C is a constant and {xz i u 3) are a particular solution 
found by integrating eqs. (|A43|) with the initial condi- 
tions 

xf(o) = 0, «f(0) = 0, 

(Xz , u%) are the solution of the homogeneous system ob- 
tained by putting D 2 = in (|A43|) and integrating up to 
R with the boundary condition 

r^O xfM~r 2 , u%(r)~2r 5 . (A45) 

For r > R eqs. (|A43|) reduce to 



The solution of this system is found as a linear combina- 
tion of the solution Xs , u f^ xt °f the homogenous system 
found by putting k 2 = J = in l|A46(l . and of a partic- 
ular solution of (|A46|) . X3 which falls to infinity 
faster than r~ 5 . Both solutions have been calculated by 
making use of MAPLE. 

The homogeneous system admits two solutions, which, 
for r — ► 00 , behave as ~ r~ 5 and ~ r 2 . Since the metric 
must be asymptotically flat, we have excluded the latter. 

Thus the general solution for r > R is 



X3 e 

U 3 C 



H 



and consequently 



^3 e 



Du 



Dwf 



w 



3ext 



where 



0) = ( 



7 



64r 3 M 7 
r 



-120r 3 M 2 2n( 
) - 45r 5 Zn(- 



r-2M 



(A47) 



(A48) 



i (A49) 



1- ■)()/• '.!//»■: 

V-2M' V-2AT 

-20rM 4 + 60r 2 M 3 - 210r 3 A/ 2 + 90r 4 Af + 8Af 5 ] 



and 



» = ( 



J 



480r 7 M' 



r ) [(1152J 2 M 9 + 700r 5 J 2 M 4 



+2100r 6 J 2 M 3 - 7350r 7 J 2 M 2 + 3150r 8 J 2 M 
+280r 4 J 2 A/ 5 + 1152Kr 3 M 10 + 64J 2 rM 8 - 320J 2 r 2 Af 7 
+ 10080.fYAf 5 r 8 + 1088XM 8 r 5 + 1664XA/V 

23520A'Af^V 7 1 fi,70n J ^ H^—^ 1 ifionnrzMSJi ' 



6720KM'r b + 168002<fMVZn(- 



2A" 



-4200r 7 J 2 M 2 ln( 
-5040iOfV/n( 
+576KM 7 r 6 ln( 

-384KM 9 r 4 ln( 



r-2M' 
r 

r - 2M 
r 



5250r*J 2 Mln( 



r-2M 



) - 134402OfVi , rc( 



r - 2Af 
r 



2Af 



) - 960 J ftTM 8 r 5 ln( 



) - 1575r y J 2 ln( 



r-2M' 
r 



r-2M 
r 



r-2M 



(A! 



The constants C and D are found by matching the inte- 
rior and the exterior solutions at r = R. 
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